Two dimensional foam rheology with viscous drag 
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We formulate and apply a continuum model that incorporates elasticity, yield stress, plasticity 
and viscous drag. It is motivated by the two-dimensional foam rheology experiments of Debregeas 
et al. [G. Debregeas, H. Tabuteau, and J.-M. di Meglio, Phys. Rev. Lett. 87, 178305 (2001)] 
and Wang et al [Y. Wang, K. Krishan, and M. Dennin, Phys. Rev. E 73, 031401 (2006)], and is 
successful in exhibiting their principal features an exponentially decaying velocity profile and strain 
localisation. Transient efi^ects are also identified. 
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While initially two dimensional (2d) foams were intro- 
duced only as a sirnple model system for numerical and 
theoretical studies [3, 13|i recent years have also seen a 
variety of rheological experiments on so-called quasi 2d 
foams, i.e. foams consisting of a single layer of bubbles 
HiQIEIHI^H]' Using bubbles trapped between two glass 
plates (Hele-Shaw cell) in a cylindrical Couette geometry 
(the foam is contained between two concentric cylinders), 
Debregeas et al. found that the flow of the foam localises 
near the inner moving wall with an exponential velocity 
profile forming shear-hands 0. While quasi static cellu- 
lar simulations 0, llDj ] , showed some agreement with the 
results, they continue to excite debate 0, especially as 
regard to the localisation of shear and deformation 0, 
which is the salient feature of the experiment. Recently, 
Wang et al. have extended shear experiments to the sim- 
pler planar geometry While their experiments using 
bubbles between a liquid pool and a glass plate showed 
the formation of shear bands with an exponential veloc- 
ity profile, a nearly linear velocity profile was obtained 
for bubble floating on the liquid (bubble raft or Bragg 
raft). This has evidenced the primordial role played by 
the method used to confine the bubbles and indicates 
that the non-uniform stress imposed by the Couette ge- 
ometry is not sufficient to explain the formation of shear 
bands with exponential decaying velocity. 

In this paper, we introduce an elementary continuum 
model for the analysis of rheological properties of a two- 
dimensional foam. It includes a viscous drag that has 
no counterpart in conventional 3d foam rheology. Our 
model is therefore closely related to the 2d viscous froth 
model llll ] which was designed to enable dynamic simu- 
lations to be undertaken with the full cellular structure 
of the foam, and included just such a viscous drag. Here 
the viscous drag will enter as a term in the continuum de- 
scription, depending on a local average of the boundary 
velocity. Experiment and theory have already addressed 
this force as it arises in the flow of bubbles in cylindri- 
cal tubes and in narrow channels j^]. It is often associ- 
ated with the name of Bretherton who showed that the 
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FIG. 1; Geometry of the case considered here in which the 
velocity v{y,t) and the displacement u{y,t) are functions of 
the y coordinate and time t. 



force varies with the two-thirds power of velocity [l^ • In 
some circumstances, a power law of one-half is suggested 
p^ . Nevertheless, as in the case of the 2d viscous froth, 
we adopt a linear form in order to keep the model and 
the analysis simple, in a search for qualitative and semi- 
quantitative understanding. In other respects, the model 
is akin to the familiar Bingham model of a substance that 
has a yield stress J4J and an internal viscosity. This, or 
one of its variants, is often invoked in the analysis of 
bulk foams. However, as in the recent work of Takeshi 
and Sekimoto jisj . we also include an elastic response, 
so that the model we propose has four key ingredients: 
elasticity up to a yield stress, plasticity, internal viscosity 
and a viscous drag force. 

While it is amenable to obvious generalisation, the 
model will be defined here for the simple planar shear 
geometry, as in 0]. Displacement u{y,t) and the veloc- 



ity v{y,t) 



9u{y,t) 

at 



are in the x direction only, as when 
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shear takes place between two parallel infinite boundaries 
in that direction (see figure QJ. This reduces the problem 
to one dimension. Strain and strain rate are reduced to 
scalars -/{y,t) = || and jiy,t) = ^. 

We will neglect inertia throughout, so that the total 
force acting on an element of fluid at y must be zero. 
Forces arise from the gradient of the shear stress a{y, t) 
and the drag force per unit area, F — —/3v, where /? 
is the mean drag coefficient. In two dimensions, stress 
has the dimension of a force divided by a length and /3 
is expressed in units of force x time per volume. The 
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required force balance is: 

da 

dy 



: (3V. 



(1) 



It remains to specify the constitutive relation for a in 
terms of 7 and 7. For simplicity, we capture the desired 
ingredients of elasticity, yield stress and plasticity with 
the following relation: 



<^ = <^Yf {i/iy) + vi- 



(2) 



Here, <jy is the yield stress and is the yield strain. 
We choose /{"//jy) = tanh(7/7y) which roughly cor- 
responds to a typical 2d static stress-strain relation for 
foams IQ . For foams 7y is of the order of unity and we 
shall set it equal to unity here. The final term in cqn. |3 
is the usual strain-rate term of the Bingham model. Note 
that for foams, the viscosity 77 depends on the strain 7. 
For low strain, the dissipation is due to the stretching of 
films and occurs at the same rate than the applied defor- 
mation. For high strain, it is mainly due to topological 
changes which leads to the disappearance and creation of 
films. This occurs at much higher rate than the applied 
deformation jl^. Nevertheless, assuming a constant vis- 
cosity is helpful in our elementary model. A very impor- 
tant restriction requires that eqn. |2] is used only when 
the strain rate 7 always has the same sign (negative in 
what follows), which is the case in the experiments we are 
referring to. In further work we will include hysteretic ef- 
fects, which are very important, but for now we accept 
this restriction. 

We can non-dimensionalise equations^and|2]by intro- 
ducing the natural length scale Lq — {ri/(3y^^ and nat- 
ural time scale Tq = rj/ay- From now, length and time 
will be expressed in units of Lq and Tq. A convenient 
representation of eqns. ^ and Q is 



where 



d'^v d J. ^ du 

dy'^ dy \ dy 



du 
'dt' 



(3) 



(4) 



The model can be solved analytically in various cases 
and limits. More generally, a numerical scheme of inte- 
gration can be used to follow the time dependence of 
the variables, as follows. We discretise y and t with 
small steps Ay and At, using lowest order expressions 
for derivatives. Given a knowledge of u in steps up to 
time t, ^ may be estimated as a backward derivative 
and eqn. 10) may be solved for v{y,t) with the imposed 
boundary conditions. Eqn. then enables us to update 
u to t + At. (In practice an Improved Euler method was 
used for the integration in time.) 

We will only consider the case in which the boundary 
at ?/ = is given a finite velocity V for all time t, while 
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FIG. 2: (a,b) Profiles of velocity v/V and magnitude of the 
strain 7 for three different times, represented by the total ap- 
plied shear T = Vt/L, shown in semi-log scale. This exempli- 
fies three regimes of exponential/linear profiles. The regimes 
also feature in (c), the variation of the magnitude of strain 
rate 7 and (d), the variation of the magnitude of stress a with 
total applied shear F, shown in log-log scale. In all the cal- 
culations shown we have chosen L = 15 and a low boundary 
velocity of V = 0.005 at t/ = 0. 



the boundary at y = L is held fixed. Correspondingly, 
u{y = 0,t) = Vt and u(y ~ L,t) — 0. For the results 
presented here, we set L = 15 and V takes various values. 
The quantity F = Vt/L may be regarded as the total 
applied shear at time t. 

The numerical results presented in Figure 2 are for 
low velocities V ^ 1 and show the existence of several 
regimes as the total applied shear is increased. 

Regime I is observed for small total applied shear, at 
which both velocity and strain profiles (Figures |21 a and 
b) are close to exponentials. Regime II is characterised 
by a linear velocity profile and a homogeneous strain. In 
regime III both velocity and strain profiles combine an 
exponential decay close to the moving boundary and a 
linear decay close to the fixed boundary. With further 
increase of total applied shear the linear tails diminish, 
leading to an asymptotic steady state (regime IV) similar 
to that for small applied shear. 

The existence of these distinct regimes is also evident 
from the plots of strain rates and stress as a function of 
total applied shear as shown in Figures [21 (c) and (d) re- 
spectively. While regime I is characterised by a strong 
localisation of both strain-rate and stress, in regime II 
(10^^ < F < 1) the strain-rate is homogeneous. The 
asymptotic steady state of regime IV is again charac- 
terised by strong localisation of the strain rate. Stress is 
saturating to its maximal magnitude 1 for all values of y. 

Figure |21 shows the velocity profiles obtained for the 
same total applied shear but for different shearing ve- 
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FIG. 3: Velocity profiles scaled by V for L — 15 (numerical 
results for □ y = 0.001, t> V = 0.3 and • F = 60). (a) 
is for total applied shear F = 0.1 and shows the succession 
of regime II, regime I and finally regime IV. (b) corresponds 
to the transition from regime II to III to IV as obtained for 
F = 1. The dashed line represents the linear solution given by 
eqn. ||SJ| and the solid line is the steady states solution given 
by eqn. Q corresponding to the exponential localisation. 



locity V. For V 1, the velocity varies linearly, cor- 
responding to the regime II. For V ^ 1, the profile ap- 
proaches the asymptotic form (regime IV). For V ~ 1, 
we can have either regime I (for F = 0.1 on figure |3| a) or 
regime III (for F = 1 on figure|21b) where we see an initial 
exponential decay followed by a linear tail (regime III). 
Figure 21 represents the different regimes encountered on 
a semi-quantitative T — V diagram. Depending on the 
shearing velocity V, several scenarios are possible before 
reaching the steady state of the regime IV. 

In order to understand these features, we return to the 
governing eqn. Q, and reduce it by various approxi- 
mations. For small time (regime I), u is small and we 
neglect the right hand side of eqn.||2J) which is approxi- 
The remaining equation 



mately — fp^ 



u = 0, 



has the elementary solution 

,sinh(L — y) 



sinh(i) 



(5) 



(6) 



Note that this solution does not vary with time, imply- 
ing that the system jumps instantaneously to the above 
velocity profile. This is indeed consistent with what is 
found in the numerical treatment, and is a consequence 
of the singular initial condition and the neglect of iner- 
tia. Provided L >> 1 in the reduced units, this solution 
is approximately an exponential over most of the range. 
The exponential profile survives until becomes large, 
and overtakes the term proportional to v. 

Neglecting the term proportional to v in eqn. |31 rather 
than that on the right hand side, and approximating the 
latter as already stated, we obtain 

d'^v d'^u 
dy^ dy^ 



(7) 



Hence u + v = a{t)y + b{t). Writing v = 
tegrating again gives u = A(y)e~* -I- a{t)y 



and in- 
+ b{t). This 
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FIG. 4: Qualitatively different velocity profiles are found 
in difi'erent regions of the F — V diagram. Regime I: expo- 
nential. Regime II linear. Regime III: combined exponen- 
tial/linear. Region IV: approach to final steady state velocity 
profile with an exponential localisation. The boundaries Fj ^jj 
and Tjj/jjj (dashed lines) are defined as in the text. The area 
to the left of the computed o data points is defined so that 
the relative error between the velocity profile and the linear 
profile characterizing Regime II (eqn. [SJ is smaller than 1%. 
A similar 1% threshold based on equation |S| has been used 
as a numerical criterion for the computation of the boundary 
between regime III and IV (□). This threshold is also used 
for the computation of the ■ data points. 



shows that the solution which develops after some time 
is linear, with a decaying transient part. The decay time 
is unity, in the units used. Applying the boundary con- 
ditions at y = and y = L, the linear variation of the 
velocity is then given by 



V 



(8) 



in excellent agreement with the simulation (see fig. ISJ- 

A further transition (regime III) takes place when the 
approximation tanh z ~ z fails, and can be replaced by 
tanh ^ 1, as the strain 7 increases beyond the yield strain 
7v'. At any given time in this regime, the second approx- 
imation replaces the first for y > y^. Thus the same 
exponential of regime I is to be expected for y < yo, 
continued by the linear solution of regime II for y > yo- 
As the time t tends to infinity (regime IV), yo tends to 
L and the solution returns to the effectively exponential 
form of eqn.ljBJ. The simulations are in excellent agree- 
ment with this profile, as shown by the solid line in figure 
01 This solution is only asymptotically reached because 
the strain in the vicinity of the fixed wall tends to zero 
due to the strong localization that appears close to the 
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moving wall. A closer analysis of this approach is possi- 
ble but will not be pursued here: suffice it to say that it 
is a slow (power law) convergence. 

The boundaries between these regions may be iden- 
tified as follows. That between regimes I and 11 may 
be found by comparing the magnitudes of the terms ne- 
glected in their respective approximations. Using the so- 
lutions given in eqns. (O and gives F//// ~ V/L, 
in agreement with the linearity in V found in numerical 
computation of the I/II boundary (see fig. ^ ). Simi- 
larly, we enter regime III when the maximum value of 
strain 7 reaches 7y, which for the linear solution occurs 
at Tjj/iji ~ 1, which is in reasonable agreement with the 
numerical data shown in figEl Putting these together, 
we see that regime II is eliminated entirely for V > L 
in dimensionless units. Reinstalling physical units this 
corresponds to a shear rate which exceeds Oy/rj. 

The steady state obtained at high applied shear F of- 
fers a very elementary candidate for the explanation of 
the phenomenon of localisation with exponential veloc- 
ity profiles in 2d foams 0, 0|. Our results allow for a 
direct comparison with the planar shear experiment on 
2d foams of ref. j^]. As eqn. ((BJ can be approximated 
(in physical units) by v/V ~ exp(— y/Lo), the velocity 
measurements provide us with a direct determination of 
Lq = (77//?)^/^. Expressed in units of bubble diameter d, 
they correspond to io ~ for the bubbles trapped be- 
tween a glass plate and a pool of liquid. Our model also 
explains why no exponential localisation is found in the 
experiments using bubble rafts. In this case, the mean 
drag coefficient /3 is expected to be very small, since there 
are no rigid plates, but rather the foam slides on under- 
lying liquid. The decay length of the exponential which 
scales like (rj/PY^^ increases up to a value of the same 
order of magnitude than the system size L and the ve- 
locity profile appears to be very close to a linear form. 
From eqn. Elwe also see, consistent with the experiments 
of 0, that in the steady state regime the scaled velocity 
profiles v/V do not depend on the shear rate V/L. 

Clearly the model can be applied more generally, for 
example to the circular Couette geometry. This suggest 
the use of polar coordinates (r, 9) which leads to an ex- 
tra term a/r in the divergence of the stress of eqn.Q. 
Although a full mathematical treatment is required to 
solve the problem in the general case, it is possible to 
use our present results, provided this extra term is much 
smaller than the viscous drag term Pv. Assuming that 
the stress is dominated by the viscous contribution 777 
during the steady state, one find that eqn. still holds 
if the distance between the two cylinders is much bigger 
than Lq. This is the case when bubbles are confined in a 



Hele-Shaw cell where the velocity profiles are found to be 
exponential with d < Lq < 2d On the contrary, for 
a bubble raft sheared between two concentric cylinders, 
the velocity profile is not exponential but rather discon- 
tinuous 0. This can be explained by a viscous drag to 
small too overcome the non uniform stress effect of the 
Couette geometry. 

We thus have seen that the exponentially decaying ve- 
locity profile in 2d foams, which is the signature of shear 
bands, is due the viscous drag generated by the bubbles 
on the confining plate. In due course more realistic forces 
(e.g. the Bretherton form) may be required, at the ex- 
pense of the extreme simplicity of what we have shown 
here. Most of the qualitative conclusions are likely to 
remain intact. 
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